function [L,U,p]=LU_factorization(A)
%[L,U,P]=LU_factorization(A)
%A: a square matrix which is nonsingular
%L: the lower-triangular matrix of A
%U: the upper-triangular matrix of A
%P: the permutation matrix records the various interchanges used

% if A is nonsigular
A_size=size(A);
L=zeros(A_size);
U=zeros(A_size);
p=eye(A_size);
A_temp=A;%save the input matrix
%check the input
if A_size(1)~=A_size(2)
    error('the matrix for LU Factorization should be square');
    return;
end
if ~det(A)
    error('the matrix for LU Factorization should be nonsigular');
    return;
end

row_temp=zeros(1,A_size(2));
pivot=0;
scale=0;
k=uint16(1);
j=k;
for j=1:A_size(2)-1
    %find the pivot, (max of this cloumn)
    [max_value,max_index]=get_abs_max_in_cloumn(A(j:A_size(1),j));
    %interchange the row
    if max_index ~= 1
        %A
        A=interchange_matrix_row(A,j,j+max_index-1);
        %p
        p=interchange_matrix_row(p,j,j+max_index-1);
        %L
        L=interchange_matrix_row(L,j,j+max_index-1);
    end
    %Gauss
    pivot=max_value;
    for k=j+1:A_size(2)
        scale=A(k,j)/pivot;
        %record in L
        L(k,j)=scale;
        %elimation
        A(k,:)=A(k,:)-scale*A(j,:);
    end
end
U=A;
L=L+eye(A_size);
%p=p
end

function [max_value,max_index]=get_abs_max_in_cloumn(A_cloumn)
%get the max absolute value in a cloumn vector
A_size=size(A_cloumn);    
if ~A_size(1)
    error('get_max_in_cloumn: input matrix is empty');
    max_value=0;
    max_index=0;
    return;
end
max_value=A_cloumn(1,1);
max_index=1;
for j=2:A_size(1)
    if abs(A_cloumn(j,1))>abs(max_value)
        max_value=A_cloumn(j,1);
        max_index=j;
    end
end
end

function A_interchange=interchange_matrix_row(A,j,k)
%interchange between two rows in a matrix
row_temp=A(j,:);
A(j,:)=A(k,:);
A(k,:)=row_temp;
A_interchange=A;
end